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Fast  Voronoi  Diagrams  and  Offsets 
on  Triangulated  Surfaces 


Ron  Kimmel  and  James  A.  Sethian 


Abstract.  We  apply  the  Fast  Marching  Method  on  triangulated  domains 
to  efficiently  compute  Voronoi  diagrams  and  offset  curves  on  triangulated 
manifolds.  The  computational  complexity  of  the  proposed  algorithm  is 
optimal,  0(M  log  M),  where  M is  the  number  of  vertices.  The  algorithm 
also  applies  to  weighted  domains  in  which  a different  cost  is  assigned  to 
each  surface  point. 


§1.  Introduction 

Voronoi  diagrams  play  important  roles  in  many  research  fields  such  as  robotic 
navigation  and  control,  image  processing,  computer  graphics,  computational 
geometry,  pattern  recognition,  and  computer  vision.  Its  Euclidean  version, 
for  which  there  is  an  efficient  implementation,  is  a building  block  in  many 
applications. 

The  Voronoi  diagram  sets  boundciries  between  a given  set  of  source  points, 
and  splits  the  domain  into  regions  such  that  each  region  corresponds  to  the 
closest  neighborhood  of  a source  point  from  the  given  set.  Let  our  domain 
be  D,  let  the  set  of  given  n points  be  {pj  € D,j  £ 0,  ..,n  — 1},  and  the 
distance  between  two  points  p,q  £ D he  d{p,q).  Then  the  Voronoi  region  G* 
corresponds  to  the  set  of  points  p £ D such  that  d{p,pi)  < d{p,pj),\/j  ^ i. 

Offsets  computation  is  often  used  in  approximation  and  singularity  theo- 
ries, and  comes  into  practice  in  computer  aided  design  (CAD)  and  numerical 
control  (NC  machines).  Given  a curve  and  its  embedding  space,  an  offset 
curve  is  defined  by  a set  of  points  with  a given  fixed  distance  from  the  original 
curve. 

There  are  some  numerical  and  topological  difficulties,  even  in  the  com- 
putation of  offsets  for  curves  in  the  2D  Euclidean  plane,  e.g.  the  formation 
of  singularities  in  the  curvature,  self  intersection  of  the  offsetting  curve,  and 
the  fact  that  an  offset  of  a polynomial  parametrized  curve  is  not  necessarily 
polynomial.  Some  of  the  numerical  difficulties  were  addressed  in  [9],  where 
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the  Osher-Sethian  level  set  method  [16,20],  which  grew  out  of  Sethian’s  ear- 
lier work  on  curve  evolution,  see  [21],  was  used  to  overcome  the  topological 
changes. 

Efficient  construction  of  distance  maps,  minimal  geodesics,  Voronoi  dia- 
grams, and  offset  curves  for  non-flat  and  weighted  domains  is  a challenging 
problem,  see  e.g.  [15,13,8,12,6,10].  The  core  of  our  approach  is  Sethian’s  Fast 
Marching  Method,  [22,19,20]  which  solves  the  Eikonal  equation  on  a rectan- 
gular orthogonal  mesh  in  0(M  log  M)  steps,  where  M is  the  total  number  of 
grid  points.  Contingent  upon  the  triangulated  upwind  and  monotonic  update 
schemes  given  by  Barth  and  Sethian  [1],  this  technique  was  extended  to  trian- 
gulated surfaces  by  Kimmel  and  Sethian  in  [11].  The  triangulated  version  of 
the  Fast  Marching  Method  has  the  same  computational  complexity,  and  solves 
the  Eikonal  equation  on  triangulated  domains  in  0{M  log  M)  steps,  where  M 
is  the  number  of  vertices.  Using  this  technique,  one  can  compute  distances  on 
curved  manifolds  with  local  weights.  For  other  applications  which  rely  on  the 
Fast  Marching  Method,  see  [14,4]. 

Here  we  apply  our  method  to  compute  Voronoi  diagrams  of  a given  set  of 
points  (or  regions),  and  to  find  offsets  from  curves  and  points  on  triangulated 
manifolds.  The  computational  complexity  of  the  proposed  algorithm  is  opti- 
mal 0{M  log  M),  its  implementation  is  simple,  and  it  also  applies  to  weighted 
domains  in  which  a different  cost  is  assigned  to  each  surface  point. 

The  key  idea  is  based  on  upwind  finite  difference  operators  as  numerically 
consistent  approximation  to  the  differential  operators  in  the  Eikonal  equation. 
Such  an  approximation  selects  the  correct  viscosity  solution.  The  upwind 
operators  allow  us  to  construct  a solution  to  the  Eikonal  equation  by  optimally 
sorting  the  updated  points  using  a heap  structure. 

The  outline  of  this  paper  is  as  follows.  The  key  for  fast  computation  of 
offsets  and  Voronoi  diagrams  is  a fast  algorithm  for  computing  the  distance. 
Hence,  we  first  comment  on  the  connection  between  the  Eikonal  equation  and 
distance  maps  on  weighted  domains.  We  refer  the  reader  to  Sethian’s  Fast 
Marching  Method  for  solving  the  Eikonal  equation  and  for  computing  dis- 
tance maps  on  orthogonal  grids,  and  to  [11]  for  details  on  our  extension  for 
computing  the  solution  on  triangulated  domains.  We  then  apply  the  method 
for  the  computation  of  fast  Voronoi  diagrams  and  offsets  on  triangulated  man- 
ifolds. 


§2.  Fast  Marching  Method  and  the  Eikonal  Equation 

We  first  explore  some  aspects  of  distance  computation  on  weighted  domains. 
In  order  to  compute  the  distance  between  two  points,  we  need  to  define  a 
measure  of  length.  A definition  of  an  arclength  allows  us  to  measure  distance 
by  integrating  the  arclength  along  a curve  connecting  two  points.  The  distance 
between  the  points  corresponds  to  the  length  of  the  shortest  curve  connecting 
them. 

Given  a 2D  weighted  flat  domain,  or  in  other  words  an  isotropic  nonhomo- 
geneous  domain,  the  distance  may  be  defined  via  the  arclength  definition.  For 
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example,  the  arclength  may  be  written  as  a function  of  the  x and  y Cartesian 
coordinates  of  the  planar  domain 

ds^  = F{x,  y)^{dx^  + dy^), 

where  T{x,  y)  : R?  —>■  is  a function  that  defines  a weight  for  each  point  in 
the  domain. 

The  distance  map  T{x,y)  from  a given  point  po  assigns  a scalar  value 
to  each  point  in  the  domain  that  corresponds  to  its  distance  from  po-  It  is 
easy  to  show,  see  e.g.  [2],  that  the  gradient  magnitude  of  the  distance  map  is 
proportional  to  the  weight  function  at  each  point 

|VT(x,7/)|  =R{x,y), 

where  |VT|  = + '^y-  This  equation  is  known  as  the  Eikonal  equation. 

The  ‘viscosity’  solution  to  the  Eikonal  equation  coupled  with  the  boundary 
condition  T{po)  = Q results  in  the  desired  distance  map. 

Our  first  goal  is  to  solve  the  Eikonal  equation.  The  key  is  to  construct 
a numerical  approximation  to  the  gradient  magnitude  that  selects  an  appro- 
priate ‘weak  solution’.  Consider  the  following  upwind  approximation  to  the 
gradient,  given  by 

(max(Z?,7T,  -D+=^T,0)^  + max(T»,7T,  -D+^T,  = Fij, 

where  for  example  D^j^T  = jg  the  standard  backwards  derivative 

approximation,  and  Ty  = T{iAx,jAy).  The  use  of  upwind  schemes  in  hyper- 
bolic equations  is  well  known,  see  for  example,  Godunov’s  paper  from  1959 
[7].  For  Hamilton- Jacobi  equations,  see  e.g.  [17,3]. 

The  solution  T can  be  systematically  constructed  in  an  upwind  fashion. 
The  upwind  difference  approximation  of  the  above  equation  means  that  infor- 
mation propagates  one  way  from  smaller  values  of  T to  larger  values.  The  Fast 
Marching  Method  exploits  this  order  of  events.  A point  gets  updated  only  by 
points  with  smaller  values.  This  ‘monotone  property’  allows  us  to  keep  a front 
of  candidate  points  that  tracks  the  flow  of  information,  ordered  in  a heap  tree 
structure  in  which  the  root  is  always  the  smallest  value.  An  update  of  an 
element  in  the  heap  tree  is  done  in  O(logM)  operations.  Thereby,  the  total 
computational  complexity  is  0(M  log  M).  We  refer  to  [22,19,20]  for  further 
details  on  the  Fast  Marching  Method. 

One  could  recognize  similarity  to  Dijkstra’s  method  [5,18]  that  computes 
minimum  costs  of  paths  on  networks.  Dijkstra  algorithm  would  obviously 
fail  to  consistently  solve  our  geometric  problems.  Actually,  any  graph-search- 
based  algorithm  induces  the  artificial  metric  imposed  by  the  graph  network, 
and  would  be  inconsistent  with  the  continuous  case,  and  thus  fail  to  converge 
as  the  graph  resolution  is  refined. 

The  Fast  Marching  Method  that  works  for  orthogonal  grids  may  be  viewed 
as  a selection  for  the  update  of  one  of  the  four  right  angle  triangles  that  share 
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the  same  vertex.  The  extension  to  triangulated  domains  is  motivated  by  this 
observation,  by  the  geometric  interpretation  of  the  update  step,  and  by  an 
additional  special  treatment  of  obtuse  angles.  We  refer  to  [11]  for  details  on 
the  extension  of  the  fast  marching  method  to  triangulated  domains.  It  is  also 
based  on  a finite  difference  approximation  to  the  Eikonal  equation,  this  time 
on  the  surface,  monotone  by  construction,  consistent,  upwind,  and  converges 
to  the  viscosity  solution. 

§3.  Offsets  and  Voronoi  Diagrams 

We  have  an  algorithm  to  compute  distances  on  triangulated  manifolds,  and 
hence  construct  offset  curves.  First,  we  solve  the  Eikonal  equation  with  speed 
.F  = 1 on  the  triangulated  surface  to  compute  the  distance  from  a source 
point  or  a region  that  defines  an  initial  curve.  We  then  find  the  equal  geodesic 
distance  curves  on  the  surface  by  interpolating  the  intersection  with  a constant 
threshold  using  a ‘marching  triangle’  procedure,  again  an  0{M)  operation. 
The  offsets  on  the  triangulated  surface,  or  the  equi-geodesic-distance  curves, 
are  shown  in  Figure  1.  The  black  curve  is  the  original  curve,  and  the  white 
curves  are  the  offsets. 

Figure  2 presents  Voronoi  diagrams  on  several  beads  and  a synthetic  head. 
We  first  compute  the  distance  from  each  of  the  initial  given  source  points 
simultaneously  using  a single  heap  structure,  and  allow  one  vertex  overlap 
between  distance  maps  form  different  sources.  The  complexity  for  the  distance 
computation  is  still  O(MlogM).  Next,  we  ‘march’  along  the  triangles,  and 
for  each  triangle  linearly  interpolate  the  intersection  curve  between  the  two 
different  distance  maps,  again  an  0{M)  operation. 

The  algorithm  complexity  remains  the  same  as  we  add  weights  to  the 
surface.  In  Figures  3 and  4 a different  cost  is  aissigned  to  each  vertex.  The 
cost,  or  weight  function,  is  texture  mapped  onto  the  triangulated  surface.  The 
weighted  offsets,  or  weighted  equal  geodesic  distance  contours  are  shown  in 
Figure  3,  while  weighted  geodesic  Voronoi  diagrams  for  several  surfaces  are 
presented  in  Figure  4.  In  both  examples,  dark  intensity  mapped  onto  the 
surface  indicates  a low  cost,  and  the  brighter  the  intensity  the  higher  the  cost. 
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Fig.  1.  OfTsets  on  four  beads  and  a Synthetic  Head. 
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Fig.  2.  Voronoi  diagrams  of  five  points  on  four  beads  and  a Synthetic  Head. 
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Fig.  3.  Weighted  offsets  on  four  beads  and  a Synthetic  Head. 


Fig.  4.  Weighted  Voronoi  diagrams  of  five  points  on  four  beads  and  a Synthetic  Head. 
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